Phase transition and annealing in quantum random energy models 
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By using a previously established exact characterization of the ground state of random potential 
systems in the thermodynamic limit, we determine the ground and first excited energy levels of 
quantum random energy models, discrete and continuous, and rigorously establish the existence of 
a first order quantum phase transition. Our analysis, corroborated by Monte Carlo simulations, 
suggests that the presence of an exponentially vanishing minimal gap at the transition does not 
prevent the success of a quantum annealing algorithm. 
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The perspective to realize a physical device represent- 
ing a quantum computer (QC) has motivated a fervent 
research activity concerning the algorithms that could 
best exploit the intrinsic quantum properties of such a 
computer as opposed to the classical ones of current com- 
puters. In particular, there has been a growing interest 
toward the possibility to use quantum annealing (QAn) 
1 as an alternative to simulated thermal annealing 
4j. A pictorial viewpoint in fact suggests that in order 
to get the ground state (GS) of a given classical Hamil- 
tonian V, the thermal fluctuations, introduced to avoid 
the system to be trapped in local minima, could be re- 
placed by quantum fluctuations able to outperform the 
former due to tunneling effects. Usually, QAn is asso- 
ciated with quantum adiabatic (QAd) algorithms 0-0] ■ 
The idea is to implement an interpolating Hamiltonian 
H(T) = V + TK, where K is an operator which does not 
commute with V. The adiabatic theorem ensures that 
for sufficiently slow changes of the parameter T the in- 
terpolating system remains in its GS so that the original 
GS of V can be recovered in the limit T — > 0. How- 
ever, for many interesting problems V the interpolating 
Hamiltonian H(T) is likely to undergo a first order phase 
transition at some value T = r c , where the energy gap 
A between the first excited state (FES) and the GS be- 
comes exponentially small in the system size N^, Q. In 
this case, a QAd decrease of T starting from some value 
r > r c , where the GS of H(T) is found with ease, re- 
quires an exponentially long time. Otherwise, the sys- 
tem evolves into a combination of several instantaneous 
eigenstates of H(T) and when r — > there is a finite 
probability to attain a state of V different form the GS. 
For V with a glassy energy landscape, "quantum is bet- 
ter" may be untrue (lp| . 

We will show that, despite the above pessimistic sce- 
nario, a QC could still be used to almost solve the prob- 
lem of finding the nearly degenerate GS of a classical 
random problem V by a QAn algorithm running in a 
time polynomial in the system size N. By this we mean 



the following. For any fixed size N, let the problem V 
represent an ensemble of random systems characterized 
by a distribution qjv(A) of the energy gap A between 
GS and FES. Unless <?jv(A) takes a Dirac-5-like contri- 
bution in A = 0, the fraction of system samples hav- 
ing A < t -1 , namely /jv(t) = J" qjv(A)dA, vanishes 
for t — > oo. By performing k annealing steps, each one 
employing a quantum computation time t, and with k 
growing no more than polynomially in N, our algorithm 
succeeds in finding the GS for a fraction 1 — /jv(t) of 
systems in the ensemble V. Note that we don't know a 
priori if the system sample considered is one of those for 
which the GS can be found exactly or not. In any case, 
however, the error in the found GS energy is 0(t ). How 
fast /at(£) approaches by increasing t is a notion which 
could be used to discriminate between simple and hard 
almost-solvable problems. Below, we discuss two exam- 
ples. In one /jv(t) < 9iv(0)f _1 , in the other /jv(i) = rjv 
constant and the problem is not almost solvable. We have 
checked these results by Monte Carlo (MC) simulations. 

We illustrate our idea in the case in which V is a gen- 
eral random energy model (REM) with discrete or con- 
tinuous distribution of the levels. For any choice of K, 
provided that K has zero diagonal elements in the rep- 
resentation in which V is diagonal, H(T) = V + TK is 
a quantum random energy model (QREM) belonging to 
the class of systems studied in [llj. For these systems, 
we have exactly characterized the GS and found a suffi- 
cient condition for the existence of a first order quantum 
phase transition Here we extend the results of 11 1 



to study both GS and FES energies of a generic QREM. 
The REM wit a Gaussian distribution is a well known toy 
model for spin glasses [l2| and the corresponding QREM 
has been first studied perturbatively and numerically in 
Q. This model provides an example of an almost solv- 
able problem. The closely related QREM with binomial 
distribution of levels turns out to be unsolvablc. 
Ground state of random potential systems. In 
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we 

have determined the exact GS of a class of Hamiltonian 
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models defined by a generic hopping operator of dimen- 
sion M and a potential taking M i.i.d. random values. 
For this class of models, in the thermodynamic limit the 
energy Eq of the GS is related to the lowest level E^ of 
the hopping operator by 



P(V) 
E -V 



dV = 



1 



E 



(or 



e < E(y (1) ), (l) 



where P(V) is the distribution of the random potential 
values V and E(V^)js the expectation of the associated 
fc-th order statistic (l3|, i.e. the expectation of the fc-th 
smallest value among the M values of V drawn according 
to P(V). Note that, whereas E^ is deterministic, Eq is 
random and ([1) is an equation for the expectation E(£o). 
However, we assume that Eq is self-averaging which jus- 
tifies the above notation. The integration interval in (|Tj) 
corresponds to the support (V m i n , V max ) of P(V). If, by 
changing a parameter below some threshold, e.g. T < T c 
in the case H(T) = V + TK, the solution of the integral 
equation in (TTJ) becomes incompatible with the constraint 
Eq < E(V(i)), then the GS of the system freezes into that 
for which V is the minimum occurred value of the given 
sample and E (T) = E(Vm) for any r < T c . A sufficient 
condition for this first order quantum phase transition to 
take place is that P(V) -t when V -> V m i D [ll|. If 
V m m = —CO) the normalization of P(V) ensures that this 
condition is always satisfied. 

Discrete and continuous QREMs. Let us consider the 
Hamiltonian H = TK + V acting on the M = 2 N spin 
states \n) = |m, n 2 , ■ ■ . , fijy) = \ni)\n 2 } ■ ■ - In^}, where 
|rij) is eigenstate of the Pauli matrix of, i = 1, . . . , TV. 
The hopping operator K is chosen as the sum of single- 
flip operators K = - Y^=i ® I2 ® • • • ® of (8 • • • ijv, of 
being the Pauli matrix acting on the i-th spin. More gen- 
eral forms for K can be considered and tackled in a simi- 
lar way. The potential V is a diagonal operator with ele- 
ments V n = (n\V\n) i.i.d. random variables drawn with 
distribution P(V). As an example of discrete QREMs 
we will consider the case in which V takes the values 
Vk = —J{N — 2k), k = 0, 1, . . . , TV, with probability mass 
function P(Vfc) = 2- Ar (^). We will refer to this model 
as the binomial QREM. As an example of continuous 
models we will assume V € (—00, +00) with probability 
density function P{V) = (TriVJ 2 )" 1 / 2 exp(-V 2 / (N J 2 )) . 
We will refer to this model as the Gaussian QREM. 

In the thermodynamic limit N — > 00, the integral in 
Eq. ([1]) can be performed exactly by the saddle-point 
method. The spectrum of TK is trivial and, in particular, 



its GS level is E, 



(0) 



-TN. Equation ([!]) thus provides 



E n 



E(V {1) ) 
-TN 



r < r c 
r>r. 



(2) 



where T c = limjv-yoo ~ E(V( 1 ))/7V. Stated in this form, 
the result applies to any QREM. In the binomial case 
E(V(i)) = — JN + 0(1), whereas in the Gaussian model 



E (^(i)) = -JNVhi2 + O(l). The critical hopping pa- 
rameter separating the frozen and paramagnetic phases 
is then r c = J and T c = J\/ln 2, respectively. 

By introducing auxiliary Hamiltonians, we can use 
again Eq. ([T]) to evaluate excited eigenvalues of H. Sup- 
pose, at first, that we are in the frozen phase TN < 
— E(V(i)). For a given realization of the random poten- 
tial V, the GS of H is \n{), where n.& is the config- 
uration of the fc-th smallest value of the potential, i.e. 
V ni < V n2 < . . . . We then introduce the Hamiltonian 
H = H — V ni |ni)(ni| whose lowest eigenvalue Eq co- 
incides with the level E\ of H . Note that H describes 
a system with random potential V having distribution 
P(V) such that E(V {1) ) = E(V (2 )). Moreover, the op- 
erators H and H have the same non-diagonal part, i.e. 



E^ 0) = E U> , so that from © we find 



p(0) 



E °-< -TN 



T < 
T > 



-E(V {2) )/N 
-E(V (2) )/JV 



(3) 



We proceed with a similar analysis in the paramagnetic 
phase. For TN S> — E(Vm), the GS of H approaches 
the GS of K, namely |0jr) = 2 _Ar / 2 E n I")- Observ- 
ing that E((0 K \H\0 K }) = -TN exactly for any T, we 
assume \0k) to represent an effective GS in all the para- 
magnetic region. We then introduce the Hamiltonian 
H = H + TN\0k)(0k\ which, we guess, has the average 
lowest eigenvalue Eq coincident with the average level 
Ei of H . Note that H describes a system with random 
potential V = V + TN2- N having distribution P(V) = 
P(V). The non-diagonal part of H is TK + TN2- N U 
where U has matrix elements U n ^ n > = 1 — 5 n , n > . It can 
be shown that any one of the TV-degenerate FESs of K, 
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A' 



1, . . . , N, is a GS of K + N2~ N U 

h(0) _ 



with eigenvalue —(N — 2) — N2 . By inserting E { 
-T{N-2 + N2- N ), and E(V" (1) ) = E{V (1} ) + TN2~ N into 
Eq. (JTJ for H, we get 



Eq = 




TN2 



-N 



T < 
T > 



-E(V(i)) 



N-2+N2-' 

-E(V ( i)) 
N-2+N2-' 



(4) 



Now we combine ([3]) and (U]) to obtain E\ . Observing 
that lirn/v_ i . 00 (E(V( 2 )) — E(Vm)) /N = 0, we conclude 
that in the thermodynamic limit 



Ei = 



E{V {2) ) T < T c 
-T(N-2) T > T c 



(5) 



Equation |T]) and hence Eq and E\ found above are 
exact up to terms O(l) [ll|. However, the O(l) errors 
obtained for Eq and E\ match when TV — > 00. It follows 
that for TV — > 00 the average gap E(A) = E\ — Eq is 
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E(V (2) ) 

2r 



E(V(i)) 
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r>r r 



(G) 



3 




r/j 



FIG. 1: Average gap E(A) in the Gaussian QREM as a 
function of V for different system sizes 8 < N < 22. The 
dashed line is the predicted value in the thermodynamic limit 
whereas for each size TV the symbols joined by straight lines 
are the averages from exact numerical diagonalization of 1000 
system samples. Inset: histograms of the frequencies of the 
gap values evaluated at T — F c (size 8 < TV < 22 from top to 
bottom, curves shifted for the sake of clarity). 

In the frozen phase, the average gap is a constant which 
amounts to E(A)/J = 2/e in the binomial case and 
E(A)/J ~ 0.6 in the Gaussian model. 

The results of Eqs. ([5]) and © have been com- 
pared with those from exact numerical diagonalizations 
of H for different values of the system size TV. In Fig. 
[T] we show the behavior of the average gap in the Gaus- 
sian QREM. In both binomial and Gaussian models, the 
data obtained for TV from 8 to 22 exhibit a systematic 
convergence toward the thermodynamic limit (|6|). 

Our numerical data confirm that, as usually assumed, 
the extensive quantities Eq and E% are self-averaging. 
For the gap we observe a different behavior. In the para- 
magnetic phase A is self averaging whereas this prop- 
erty is lost in the frozen phase. In the inset of Fig. Q] 
we show the frequencies of the gap values at the crit- 
ical point for the Gaussian QREM. It is evident that 
limjv^oo var(A)/E(A) 2 ^ 0. The same conclusion is 
reached for both binomial and Gaussian QREMs in the 
whole range < T < T c . 

The gap distributions for the binomial and the Gaus- 
sian QREMs have well distinct shapes in the frozen 
phase. In the discrete model, when T is decreased below 
r c the gap values concentrate around A = 0, 2 J, 4 J, .... 
Including only terms not vanishing in the limit TV — > 
oo, at r = the gap distribution becomes gjv(A) = 
rjy(5(A) + (1 — rjy)5(A — 2 J), where can be calculated 
analytically and for TV large approaches 1 — 1/e. In the 
continuous model, the gap values are evenly distributed 
in the range < A < 2r c for T < T c . The distribution 
smoothly deforms by decreasing T and at T = becomes 



gjv(A) ~ gjv(0) exp[— gjv(0)A], where gjv(0) quickly sat- 
urates for TV —> oo to a constant value near 1/(0. 6 J). 

Minimal gap and QAd algorithms. To complete the 
study of the gap, it is useful to define two other quanti- 
ties [8j. For any chosen system sample of size TV, A(T) 
has a minimum at T = r m j n , see Fig. [2] Let us define 
A m in = minr>oA(r). Theoretical arguments, checked 
numerically for both Gaussian and binomial QREMs, 
show that i) A m j n and r m ; n are self-averaging quan- 
tities, it) E(A min ) ~ 2T C N2~ N / 2 for TV large and Hi) 
limjv^oo E(r min ) = T c with var(r min ) ~ (T c /N) 2 for TV 
large. In fact, the shape of A(r) for a chosen system 
sample can be accurately reproduced by Ritz method 
writing the lowest eigenstates of H(T), \E n (T)) with 
n = 0,1, as a superposition of the lowest eigenstates 
of K, \0 K ) and |1^), and of the lowest eigenstates of 
V, \n{) with I = 1,2, ... , together with the associated TV 
first neighbors \n^) = af\rii). The crude approximation 
\E n (T)} = a„(r)|0if) + b n (T)\ni) provides the analytical 
estimates r min = -V n JN and A min = -2V ni 2~ N ' 2 . 
Properties i), ii) and Hi) follow straightaway from these 
expressions observing that V ni = V(i) and var(V(i)) ap- 
proaches a constant close to T 2 for N — > oo. 

Note that for r < r c , where A is not self-averaging, 
E(A) loses memory of the sample dependent minimum of 
the gap, see Fig. [TJ On the other hand, the self-averaging 
property in force for T > r c implies that E(A) ~ A 
changes smoothly, at finite N, between the frozen and 
the paramagnetic regimes for T G [r m i n ,r m i n + r c /iV]. 
For later use we also note that for |T — r min | < T c /N the 
superposition of the GS \E (T)) with both states \0k) 
and \ni) is much larger than (0jf|ni) = 2~ N / 2 . In the 
crude approximation, we have ao(r) ~ 1 for T > r m j n 
and oo(r) ~ l/[2 + (N(T-T mia )/(2V ni 2' N / 2 )) 2 }- 1 / 2 for 
r < r m i n . Identical expressions hold for bo(T) in the 
reversed regions T < r m ; n and T > T m - m , respectively. 
In the Ritz method with a larger superposition of states, 
\Eq(T)) has for |T — T m \ n \ < T c /N a substantial super- 
position also with a few of other V eigenstates, namely 
\ n 2): 1^3}; ■ • ■ j an d their associated N first neighbors. 

Points i) and ii) imply that a QAd algorithm to find 
V ni would require a computational time t > A~J n ex- 
ponentially long with TV for almost any system sample. 
In principle, one can hope to elude this limitation by us- 
ing an adiabatic path in which the transit through a gap 
minimum is avoided [14j . However, it is not clear how to 
realize this path for V of the class considered here. 

Imaginary time non adiabatic QAn algorithm. In al- 
ternative, consider the following non adiabatic strategy. 
Chosen a system sample of size TV, let it to evolve for 
an imaginary time s € [0, kt] with the time dependent 
Hamiltonian H(T(s)), where T(s) = J2j=o ^jXyt.jt+t] («), 
X[a.b] being the characteristic function of the interval 
[a, 6]. For simplicity, we set Tj = Tq — jST, j = 1, . . . , k, 
with ~ 0. We also choose Tq ^> r m ; n so that we 
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can initialize the system state to the approximate GS of 
H(T ), i.e. |V(0)) = |0jc>- If ^ ~ 2r c /A, even if V ni is 
unknown, with probability 1 for N — > oo, one the Tj val- 
ues, let say r 3 -„ , occurs at time s = j*t in the gap valley 
of the system sample [r min — T c /N, T min + T c /N}. 

At the generic time s = jt, when T(s) suddenly drops 
from to Tj, the system state \i/>(Jt)}, which has par- 
tially relaxed to the GS of H(Tj-±), docs not change. 
Nevertheless, rewritten in the basis of the eigenstates of 
H(Tj), it is a non trivial superposition which, in the next 
interval [jt,jt + t], partially relaxes to \Eo(Tj)). In a for- 
mula mt + t)) =J2 n (E n (T j )mt))e- E ^ )t \E n (T j )). 

Untill j < j* we have | J E (r j )) ~ |0jf> and \ip(Jt + 1)) 
remains close to the GS of K. At j = j*, \E n (Tj t )) has 
a substantial superposition with both \0k) and |ni). It 
follows that the probability to find the system at time 
j,t + fin the GS of V, |(ni|V>(j*i + t))\ 2 /\\ip(j,t + t)\\\ 
is much larger than the bare value |(ni|0_fs-)| 2 = 2~ N . 
As discussed before, transitions into a few of other V 
eigenstates, in particular ln.2) , \n^), • ■ • , take place with 
probabilities of similar order. The evolution steps j > , 
during which \Eo(Tj)} ~ |ni), filter out all these excited 
eigenstates of V, provided that A(0) = V n2 — V ni > i -1 . 
Since A(0) is non self-averaging, it depends on the system 
sample if this condition does apply or not for a chosen t. 
In the binomial case, the fraction of system samples of 
size N having A(0) < i -1 is /jv(f) = r\/v, whereas in the 
Gaussian model /i\r(t) ^ gjv(0)/i. We conclude that for 
N large the fraction of Gaussian samples for which the 
QAn algorithm does not succeed vanishes as £ _1 , whereas 
in the binomial case this fraction is 1 — 1/e. 

We have verified the above scenario by means of MC 
simulations. The matrix elements of the imaginary time 
evolution operator Texp(— J" fct H(T(s))ds) obtained by 
our MC algorithm [THj] are exact except for the statisti- 
cal error caused by the use of a necessarily finite number 
of random walks in the configuration space. Figure [2] 
shows the results for one Gaussian QREM sample of size 
N = 20. At the end of each time interval of length t, the 
hopping parameter is decreased stepwise from the value 
r = J with ST = 0.1J. For every step we compare the 
exact GS energy E (T), obtained by numerical diagonal- 
ization, with E$* C (T), the GS energy estimated by the 
MC executing a set of T random walks. In Fig. [2] we 
also show the gap A(r) of the system sample under ex- 
amination. This gap has a minimum A m j n ~ 0.018J at 
r m i n ~ 0.7 J, whereas its value at T = is A(0) ~ 0.06J. 
For t = 8/J, i.e. A(0) < i" 1 , the value of £^ /c (r) ob- 
tained for r — > may occasionally provide the correct 
result £7o(0) = V ni , see data for T = 2 16 , but erratic re- 
sults are produced for different values of T. On the other 
hand, whenever A(0) > t^ 1 a systematic convergence 
E^ C (T) -> E Q (0) is observed for r -> by increasing 
T. Note that for t = 32/ J the error E^ C {T) - E {T) 
is large near r m ; n where A(T) < t —1 . For t = 64/ J, i.e. 
A m in > a QAd condition is achieved and Eq IC (T) is 
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FIG. 2: Relative error in evaluating the GS energy by MC 
simulations (symbols joined by straight lines, left scale) in a 
Gaussian QREM sample with iV = 20 spins. The hopping 
parameter T is decreased stepwise from the value F — J and 
at each step Eq IC (T) is evaluated by executing T random 
walks in the configuration space for a time t — 8/J (unfilled 
symbols, T = 2 14 , . . . , 2 20 ), t = 32/J (filled symbols, T = 
2 14 , . . . , 2 20 ) or t = 64/J (crosses, T = 2 20 ). The gap of the 
system sample A(F) is also shown (solid line, right scale) in 
comparison with the chosen t -1 values (dashed arrows). 



close to Eq(T) for any value of T. 

In the MC implementation of our QAn algorithm the 
total computation time is proportional to ktT. By in- 
creasing N, provided A(0) > < _1 , £0(0) is found if the 
number of annealing steps scales as k ~ N and the num- 
ber of random walks as T ~ 2 N . Thus the total computa- 
tion time increases exponentially with the system size TV. 
The conclusion would be drastically different by using a 
QC, say with a number of qubits > N. In that case, we 
wouldn't be in need of a classical probabilistic represen- 
tation of the matrix elements of Texp(— H(T(s))ds) 
and the total computation time would drop to kt. 
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